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Abstract 

In this paper we consider a system of three fractional differential equations describ- 
ing a nonlinear reaction. Our analysis includes both analytical technique and nu- 
merical simulation. This allows us to control the efficiency of the numerical method. 
We combine the numerical approximation of fractional integral with finding zeros 
of a function of one variable. The variable gives a desired solution. The solution has 
a peak. Its position in time depends on the order of fractional derivative. The two 
other solutions of this system have a completely monotonic character. 
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1 Introduction 



The dynamical systems with fractional dynamics are adequate for the de- 
scription of various physical phenomena in viscoelasticity, diffusion processes, 
dielectric relaxation, etc [T1I21I3B] . However, now the analysis of nonlinear dif- 
ferential equations with fractional derivatives is in early development. The 
point is that on the one hand the derivation of such equations from some ini- 
tial principles is in the making. On the other hand the numerical simulation of 
fractional differential equations are enough expensive |5]. In the present state 
of affairs it may follow a formal way by changing the derivatives of integer 
order in the well-known nonlinear equations on the derivative of an arbitrary 
order. This means a search to the touch together with overcoming the difficul- 
ties of numerical analysis. So that was at the beginning of fractional calculus. 
But there exists an alternative way. The other approach may be based on the 
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recent success in the study of linear fractional differential equations. In partic- 
ular, the linear fractional equations have been well investigated analytically in 
detail p[7|5] . Moreover, the fractional oscillator and the fractional relaxation 
have an clear physical interpretation PITU] . Some nonlinear models turn out to 
be very close to linear ones in a sense. This permits one to combine analytical 
computations with simple numerical simulations. The aim of this work is to 
present one of such nonlinear systems under fractional dynamics. 

The plan of the paper is as follows. Section |2] is devoted to a derivation of 
the nonlinear reaction obeying a stochastic time arrow. The arrow leads to a 
system of fractional differential equations that describe the process evolution. 
One of the equations can be resolved by analytical methods. The sum of the 
two others has also an clear analytical form. The feature analysis is represented 
in Section |3l Nevertheless, the system of the equations upon the whole is 
nonlinear. Therefore, the next step is directed to numerical simulations of those 
equations what cannot be resolved analytically. We rewrite them in the integral 
form and apply a recursive procedure of numerical fractional integration. It 
should be observed here that the results obtained in Section [3] proves to be 
useful for Section HI We sum up and discuss the outcome in Section [51 



2 Theoretical model 

The nonlinear processes play an important role in the kinetics of chemical 
reactions [TT]. The typical example is a reaction in the gas that consists of 
atoms of two grades {A and B). The reaction A + B = AB creates two-atom 
molecules. In the same time it is not impossible an inverse process, dissociation 
AB = A + B. The dynamical system can be described by differential equations 
in the nonlinear form. 

One of similar nonlinear reactions is mentioned in [12] . Its temporal evolution 
is written as 



with the initial condition N^/N = yi{0) = 1, N2/N = 7/2(0) = 0, N^/N = 
2/3(0) = 0, where N is the common number of atoms, and A^^i, N2, the 
number of molecules of three types arising from the chemical reaction of the 
atoms. The equation for yi{t) is linear and is independent from y2{t) and 
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Usit). Using the stochastic time arrow described in [9], the equation is easily 
generahzed in the fractional form. 

The stochastic time arrow characterizes an interaction of a physical system 
with environmental degrees of freedom. Then the time variable becomes a sum 
of random temporal non-negative intervals Tj. They are independent identi- 
cally distributed variables belonging to a Levy-stable distribution. Their sum 
gives asymptotically a stable distribution with the index < a < 1. Following 
the arguments of piTS] , a new time clock is defined as a continuous limit of 
the discrete counting process Nt = max{n G N : X^iLi ^« — 0) where N is 
the set of natural numbers. The time clock is the hitting time process S{t). 
Its basic properties has been represented in [TSfTl] . The probability density of 
the process S{t) is written in the form 

/(t, r) = ^ y e"*-™" m"-^ du , (2) 

Br 



where Br denotes the Bromwich path. This probability density has a clear 
physical sense. It defines the probability to be at the internal time r on the 
real time t (see details in [9]). 

Let us make use of the general kinetic equation 

^ = WpW. (3) 



where p{t) defines the transition probability (or the concentration of a sub- 
stance) in a system, W is the transition operator between the system states. 
In a general case the operator W may be even nonlinear. Provided that it 
is independent on time explicitly. Now we determine a new process with the 
probability 

oo 

Pa{t) = J p''{t,T)p{T)dT. (4) 




It is that the function p^{t,T) characterizes the stochastic time arrow. The 
Laplace transform of the probabilities Pa{t) and p{t) gives the relation Pa{s) = 
s°'~^p{s") with < a < 1, and the Laplace image of x{t) is written as 

oo 

x{s) = J e''^x{t)dt. 
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By the simple algebraic computations we have 

After the Laplace inverse transform the fractional representation of Eq.Q 
takes the form 

1 I- 

Pait) = p(0) + —- ldT{t- r)"-i Wp„(t) . (5) 

Although Eq.(l5]) is integral, it can be expressed in terms of fractional derivative 
P)13fl4] . This approach permits one to generalize Eqs.([l]) in that way too. For 
more details, see the next section. 



3 Analytical calculations 

Consider the following nonlinear system described by fractional differential 
equations [15j 

yi-y'L (6) 



=7/1 — ?■ 

dt"" 

d'^ys _ 2 



dt 



a 



2/2 



under the initial condition 

?/i(0) = l, 1/2(0) = 0, 1/3(0) = 0, 

where < a < 1. Here the fractional operator d'^/dt'^ is supposed in the sense 
of Caputo 



d"x{t) 1 f 



I : — - — dr, n — \ < a < n, 

J (t-r)"+i-" 
^ ' 



dt'' T{n -a) J {t - r)"+ 

where x^"'\t) means the n-derivative of x{t). To sum up Eqs.([6]), we get 
d"{yi + y2 + y3) 



dt' 



2/1 + 2/2 + 2/3 = 1 
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The system analysis shows above all that the first equation for yi is linear. 
It may be integrated analytically, namely yi{t) = Ea{—t°'). The variable yi is 
expressed in term of the one-parameter Mittag-Leffler function Ea{z). Thus, 
the value yi becomes vanishingly small with >-oo asi~"/r(l — a). On the 
other hand, from the second and third equations it follows that the sum of y2 
and yz gives 

This means that the definition of yi allows one to establish the functional 
dependence between y2 and y^. Therefore, the sum takes the form 

y2 + y3 = l-£^a(-n- 

This result tends to one as i — >■ oo. 
We find it useful to define the identity 



where Ea^p[z) = ^'^=oz'^/r{an + /3) is the two-parameter Mittag-Leffler func- 
tion, and Ea{z) is a particular case of Eot^p{z), namely Ea,^i{z). The Laplace 
image of t"-Ea,i+a(~^") is written as 

/ e""* i+a(-t") dt=- , 

^ ' 

while 

oo a—1 

^ ^ 



The sum of these images reduces to 1/s. 

In the upshot it still remains to find the variable y2- The value satisfies a 
nonlinear fractional differential equation with the right-hand term. In other 
words this equation is inhomogeneous: 

Its form is such that for the final conclusion it is necessary to provide a numer- 
ical simulation. Nevertheless, our analytical analysis will not be full without 
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Fig. 1. Nonlinear reaction with ordinary time evolution ([T|) . 

asymptotic estimations as applied to ?/2- If a is equal to one, we can neglect 
the right-hand term because of its exponential decay. The value 1/2 (^) behaves 
algebraic for t ^ 1. If the index a is otherwise (0 < a < 1), the right-hand 
term of Eq.Q is the one-parameter Mittag-Leffler function having an alge- 
braic asymptotic form. What term in Eq. ([7]) gives a significant contribution in 
a given case? The feature will be established by a numerical simulation more 
clearly. The variable ?/2(^) becomes vanishingly small, when t tends to infinity. 

Moreover, the Adomain's composition method assumes a series solution for 
[15] . The first three terms of this series are given by 

r t^a / r(2a + l)\ t^'" 

^'^^^ ~ r(l + a) ~ r(l + 2a) ^ 1^ ~ Y\a + 1) ) r(3a + 1) ' ' ' ' 



They indicate the presence of a maximum. However, the determination of 
its properties (position and value) is not enough reliable for the composition 
method. The point is that the Taylor series converges only for < t < 1, but 
the maximum is the case for t > 1. Thus the estimation ([HD for the maximum 
will not be exact. This can be easy to verify numerically for a = 1, using the 
usual methods (for example, the Runge-Kutta technique or something like 
that). In this case, the system ([6]) consists of first-order differential equations. 
Their numerical solutions are represented in Fig. [1]. 



6 



4 Numerical approximation 



To solve the nonlinear fractional differential equation, the principal approach 
is to use a numerical approximation of fractional integral. Its advantage is due 
to the fact that the kernel of this integral is regular for a > 1, and weekly 
singular for < a < 1, hence integrable (at least in the improper sense) for any 
index a > 0. In contrast, the kernel of fractional derivative is not integrable, 
and it requires special regularization methods restricting possibilities of the 
numerical simulation |17] . 

Let {tn = nh : n = 0,1,2, ... , N} be a partition of the interval [0 T] into a 
grid, where h = T/N. According to [T7j, the fractional integral of the value y 
is computed as 



where the quadrature weights are derived from a product trapezoidal rule to 
the following: 



(1 + a)N'^ - N'^+^ + {N- , if n = 0; 
1 , iin = N; 



(A^-n + l)"+i -2{N -n) 
[ +(A^ - n - 1)"+^ 



a+l 



+ 



if < n < N -1. 



The error of this algorithm is of order 0{h'^). The accuracy can be improved 
by the Richardson extrapolation procedure [T7] . 

As an example to this algorithm, we suggest to consider the fractional one- 
half-order integral of the exponential function. The integral has an analytical 
representation: 



1 'r 

J'/' exp(-t) = y (t - r)-V2 e-^dr = 

= -— j e- sin TdT = e- *erfi( v^) , (10) 

where erfi(x) = ^ Jq exp{y'^) dy is the imaginary error function. The tran- 
scendental function is expanded into a Taylor series about x: 

2 oo ^2m+l 

erfi(a;) = —= ^ 



TT (2m + l)m! ' 
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Fig. 2. Fractional one-half-order integral of the exponential function, J^/^exp(— t), 
obtained from Eq. ijlOp (line) and by means of Eq.([9|) (circles). 

The expansion allows one to compare J^^^ exp(— t) in the analytical form (ITO!) 
with the numerical simulation result of the fractional integral by means of the 
above-mentioned algorithm ([9]) (see Fig. [2]). 

Next we transform Eq.(I7]) into the integral form: 



r 



+ J''yi{t) = ^2(0) + y2{t) + ry'^{t) = 1 - E^{-e) . 



It is useful to remark the following property: < y2{t) + J'^y^it) < 1. As 
J'^y'^it) > 0, it is easy to find < y2{t) < 1. The upper limit of y2{t) is too 
strong. As will be shown below, y2{t) is less than one. The positivity property 
of 1/2 (i) at any time is very important. The same exists also for yi{t) and 1/3 (t). 
Both functions yiit) and ysit) are completely monotonic, but y2{t) is not. The 
latter has a maximum. 

Applying the approximation ([9]), Eq. (fTT]) is written as 
r(2 + a) 



The expression takes the form of an ordinary quadratic equation: 
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Fig. 3. Nonlinear reaction with time evolution under fractional derivative with the 
index a = 0.8 . 

l,a ua N—1 

r(2 + a) T{2 + a) 

-l + E„(-iV"/i") = 0. (12) 

Then the recursive procedure is obtained in an elementary manner. The nu- 
merical scheme permits one to calculate y2{nh) from y2{{n — l)h) to solve the 
quadratic equation. Each subsequent value of 7/2 depends on previous ones, 
because the fractional integration exhibits a long-term memory loss. For com- 
puting the one-parameter Mittag-Leffler function Ea{—t") we use a numerical 
scheme listed in Algorithm 4 [T7j that is reproduced with minor corrections 
from [18] . The quadratic equation ( fT2|) has two solutions. One of them is pos- 
itive, and another is negative. We select the non-negative solution in the force 
of the above analysis to the analytical properties of 1/2 (^)- The temporal evo- 
lution of this system for a = 0.8 is shown in Fig. [3l Formally, the system 
behavior for < a < 1 is like to the system one for a = 1: the variable yi{t) 
monotonically tends to zero, the variable 2/3 (t) to one, and the variable 2/2 (^) 
possesses one maximum after which it decays to zero. But there are differences 
that we discuss in the next section. 

This approach to the numerical simulation of the system ([6]) has the following 
advantages. Firstly, the numerical errors are of the order 0{h?). If we use 
Algorithm 3 from [T7j, they would be of the order 0(/i^+°), i. e. they would be 
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Fig. 4. Peak amplitude of the variable for 0.1 < a < 1 in the nonlinear reaction 
with fractional dynamics . 

slightly more. Moreover, Algorithm 3 is not less expensive than the procedure 
([9]). A reduction in the amount of computational work can be achieved by 
using a graded mesh, thereby taking the 0(A^^) method of ([9]) and making it 
0(A^log A^). This will be the same as in Algorithm 3. Thirdly, the Richardson 
extrapolation procedure, improving the numerical simulation accuracy, does 
not have been proved still for Algorithm 3. Hence, its application to Algorithm 
3 is only based on assumptions of its expected correctness [17]. As for the 
procedure ([9]), the use of the Richardson extrapolation is exactly true for any 
a > 0. It is also worth noticing that Algorithm 3 itself cannot yield a numerical 
solution of nonlinear fractional equations. The algorithm should be added by 
calculating zeros of the corresponding function of one variable. It is that the 
variable is a required solution. However, this way to numerical simulations 
of nonlinear fractional equations is not considered in the review [17]. In our 
opinion the approach has good perspectives. 



5 Results and Discussion 

The variables y2 and 7/3 describe a relaxation of substances in some phys- 
ical system more than anything yet. And what is why. One of its possible 
interpretations lies in the representation of the process as a reaction. Some 
substance promotes an initiation of two other interacting substances. One of 
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Fig. 5. Time position of the ?/2's maximum in the dependence of the index a for the 
nonhnear reaction with fractional dynamics . 

them is metastable and transforms into another. But this reaction does not 
occur at once. It evolutes in time. Therefore the change of the variables yi, y2 
and 1/3 manifests something hke a relaxation. 

Let us choose the index a as a governing parameter. This means that we will 
change it from zero to one for establishing how it influences on the character of 
the system The key characteristics of this system is a temporal position of 
the 2/2 's maximum and its value. The results can be only obtained by numerical 
simulations. They are represented in Fig. H] and Fig. [51 The first observation 
shows that the model ([6]) with the index a = 1 is well agreed with the models 
for a 1. Otherwise speaking, all they together form a smooth line in the 
interval < a < 1. Figure Ohas been shown for 0.4 < a < 1 to be easy to 
recognize a minimum in this dependence on a. Really, the minimum is the 
case for a ~ 0.8. This indicates that at first the temporal position comes to 
the abscissa axis, when the index a changes from one to ~ 0.8. However, the 
next decay evolution of the index leads to the motion of the 1/2's maximum 
outside of the abscissa axis. The peak amplitude has a completely monotonic 
character depending on a (see Fig. H]). 

By the concrete example we have shown that the numerical procedure of 
fractional integration can provide a convergent method for the solution of 
nonlinear systems of fractional differential equations. The approximation series 
solution (Adomian decomposition) is not so successful because of the too slow 
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convergence of the corresponding series expansion for such a type of equations. 
Therefore, in this case the decomposition plays only a secondary role. Although 
the approximation series solution can be obtained to any desirable number of 
terms, this makes the method expensive and ineffective. Probably, it would 
be better to add the algorithm by an asymptotic expansion. However, it is 
difficult to suggest any universal method. Thus, the numerical simulations turn 
out to be especially useful for the analysis of nonlinear fractional differential 
equations. 
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